Theoretical Population Biology

The Yampolsky-Stoltzfus (YS) model provided the first formal demonstration that a bias in the introduction of new mutations can shape the direction of adaptive evolution, even when mutation rates are small relative to selection coefficients. Despite its foundational role, the model has been characterized almost entirely by simulation and closed-form solutions for its central quantity, the probability PA that a mutationally favored but less-fit genotype fixes ahead of a rarer, fitter alternative, have been available only in the limiting regimes of vanishing mutation supply and infinite population size. Here we derive closed-form approximations for PA across the intermediate regime of finite population size and appreciable mutation supply that lies between these two limits. Under strong selection between the competing mutants, we combine an establishment-probability argument with clonal interference theory to obtain an explicit expression for PA that accounts for both the ongoing input of the mutationally favored mutation and the frequency-dependent establishment probability of its competitor. Under weak selection between competing mutants, we derive a complementary expression based on the fixation of a mixed population under Kimura's classical formula. Both expressions agree closely with direct simulation of the underlying Wright-Fisher process, recover the known limiting cases, and remain consistent with one another over a substantial intermediate range of parameters. These results are expected to provide a basis for assessing how the strength of mutation bias affects adaptation in empirical data.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Many pathogenic and microbial species undergo complex life cycles wherein they experience drastic bottlenecks during invasion into a host and rapid growth within hosts. The effects of such a complex population history on patterns of genetic variation is not currently understood. We model a pathogen population as a metapopulation where each deme represents a host and incorporate transmission dynamics as well as the effects of transmission bottlenecks. We employ a coalescent framework, obtaining analytical expressions for pairwise times to coalescence and genetic differentiation within and between hosts. We find that recurrent bottlenecks rescale the coalescent process within hosts to the size of the bottleneck, reducing within-host variation. In addition, the pairwise distribution of time to coalescence under our model deviates substantially from the Wright-Fisher process when the number of hosts (d) is of similar order to the size of the bottleneck (k), i.e., k [~] d and when m <<1/(2k) +1/(2d), where m is the migration rate. Otherwise, the pairwise process under our model converges in distribution to the Wright-Fisher process. Coalescent simulations show that under the same conditions the site frequency spectrum (SFS) under our model deviates from that under the Kingman coalescent and is drastically skewed towards singletons. Our results suggest that the inference of selection or demographic history from pathogen population genetics data may be systematically biased by the mechanics of pathogen reproduction, indicating a need for further theoretical work that can inform inference.

Forward-in-time population genetic simulations are widely used in evolutionary analyses, but simulating large populations and long genomic regions remains computationally demanding. To reduce this cost, parameter rescaling is widely employed, in which the original evolutionary process is approximated by one with a smaller population size and fewer generations. Recently, several studies using the SLiM simulator have raised concerns about the accuracy of this rescaling approach. In this study, we show that many of the biases reported in these studies can be mitigated by using a different simulation algorithm. These results reveal that the accuracy of parameter rescaling depends on how well the simulation algorithm preserves diffusion-limit properties under rescaling.

Amongst measures of linkage disequilibrium, the correlation measure r2, due originally to Sewall Wright, has several advantages. Although basically a haploid statistic, it can readily be estimated using covariance and correlation from diploid data. Its population expectation can also be found from probability methods. We show that the probability of linked identity by descent, LIBD, directly estimates the expectation for r2 for Wright-Fisher models. The simplest such model for two loci is the haploid model, where each haplotype is chosen independently from haplotypes of the previous generation. This gives an equilibrium between recombination (recombination frequency c) and drift (population size N) where the expected value of r2 is approximately 1/(1+4Nc). For diploid models, the expectation is similar for low values of c, but contains a term in c2 and rises to double the haploid value for the limiting value of c = 0.5. An increase in r2 is also expected for the model of monogamous mating, as opposed to random mating for each offspring. The expectation for r2 is unaffected by the use of diploid rather than haploid data, an exception being the case of monogamous mating. Computer simulation shows close agreement of r2 with probability expectations, with minor discrepancies observed for unlinked loci. Using the probabilistic definition, r2 can be intuitively understood as a parameter shaping the effective sample size in genome-wide association studies.

Evolutionary rescue occurs when genetic adaptation prevents extinction of a declining population. Mutation and recombination generate adaptive variation, but when their rates at adaptive loci are too low, rescue may require evolution toward higher rates of genetic change. Here, we derive and compare analytical approximations for the probability of rescue through two alternative routes: mutation-rate and recombination-rate evolution. Our model comprises two resident genotypes from which the same adaptive haplotype can arise along either route, with rescue requiring the appearance of an unlinked modifier allele that increases mutation at, or recombination between, two fitness-determining loci. We show that resident genotype composition determines the relative success of the two routes. Mutation-rate evolution yields a higher rescue probability when one resident genotype is sufficiently more common because mutation can act on either genotype whereas recombination requires both. Otherwise, the two routes yield similar rescue probabilities. Even when either route makes rescue appreciably likely, however, the modifier alleles facilitating rescue are often lost before rescue is complete. This loss is driven by multiple factors: segregation can separate the modifier from the adaptive haplotype; the fitness cost of modifier expression raises extinction risk for modifier-bearing lineages; and continued modifier-induced mutation disrupts the adaptive haplotype after it forms. The latter effect can make recombination-rate modifiers more likely to persist even when mutation-rate evolution yields a higher rescue probability. Overall, these results show that both mutation- and recombination-rate evolution can induce rescue, but successful rescue need not entail persistence of the modifier allele that facilitated it.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Haldane's principle relates equilibrium mutation load to deleterious mutational input rather than to selection strength. We study this balance in a homeostatically regulated, damage-structured population where mutation occurs during reproduction, repair restores individuals toward lower-damage states, and nonlinear demographic feedback maintains constant population size. An exact demographic identity shows that load at every time equals excess reproductive effort weighted by mean relative fecundity. Without repair, the equilibrium branch with positive pristine-class frequency is exactly solvable, and along it the load is exactly the mutation probability times the realized per-capita reproductive turnover rate; when the terminal cap is removed, this is precisely the realized deleterious input, for every mutation probability below one. For rare mutation it reduces to the mutation probability times the pristine turnover rate, with the local fitness contrast canceling at leading order. With a finite terminal class the branch instead reaches a closure-dependent boundary collision. Repair modifies the rare-mutation balance by partitioning first-order removal from the first damaged class between selection and restoration, multiplying the leading load by the selective share of that removal. An exact binary solution extends these results to finite mutation pressure: the unrepaired interior equilibrium collides nonhyperbolically with the fully damaged boundary and converges algebraically at the threshold, whereas any positive repair rate restores a unique interior equilibrium for positive mutation; at the no-repair threshold, its displacement from the boundary scales as the square root of the repair rate as repair vanishes. These results separate the familiar competing-transition structure of mutation--selection balance from the specifically homeostatic effects of turnover-dependent mutation input and nonlinear demographic feedback.

Regular environmental cycles such as day and night or the lunar cycle are ubiquitous and drive many biological processes. Many organisms have evolved endogenous time-keeping mechanisms to anticipate these cycles to accurately time their behavior or synchronize physiological processes. The evolution of such synchronization to environmental cycles and their variability is not very well understood. As an example, we use the marine insect Clunio (Diptera, Chironomidae), which synchronize adult emergence and reproduction with the lunar cycle, and show different emergence patterns across populations. We use a model constructed with replicator-mutator dynamics to provide a novel outlook at the evolution of these emergence patterns. In particular, we investigate the effect of two traits of individual phenotypic variability of emergence, the spread and shape, on population-wide emergence patterns. We find that the semilunar and lunar emergence patterns are alternatively stable and occur when the phenotypic distribution has a narrow and is roughly flat. Arrhythmic emergence occurs when the phenotypic distribution is wide and nearly flat. The semi-lunar pattern can only be observed when the parameter of rhythmicity is changed. With our model, we introduce the concepts of frequency-dependent selection and individual phenotypic variability into the theory of chronobiology, which gives us insights into the evolution of rhythmicity in these organisms.

Demographic inference based on sequentially Markovian coalescent (SMC) model is widely used for disentangling the evolutionary history of species and populations. However, the genetic effects of self-fertilization, including the expected reduction of population-level genetic diversity and hence effective population size (Ne), are often ignored in demographic inference, despite the fact that selfing is a common reproductive mode in both invertebrates and plants. To investigate how self-fertilization could bias demographic inference, we apply eSMC2, a selfing-aware demographic inference method, to Hawaiian populations of Caenorhabditis elegans, a predominantly self-fertilizing nematode. Using an expanded dataset that contains more Hawaiian strains than in previous studies, we recovered previously reported genetic groups, including one largely consisting of strains from the Big Island. Population substructure is also detected among these Big Island subpopulations, highlighting the importance of identifying subgroups to avoid misinterpretations of demography. For different Big Island subpopulations, we estimate high selfing rates compatible with previous studies and infer consistent population declines. Population sizes inferred using eSMC2 are smaller than those inferred using PSMC', which does not account for selfing. However, rescaling PSMC' parameters to account for self-fertilization produces a similar demographic trajectory as eSMC2 inference. Our study thus demonstrates how one can account for effects of selfing when using SMC-based methods. It also indicates that the Hawaiian group shows population declines similar to those previously found outside Hawaii.

Oncogene amplification on circular extrachromosomal DNA (ecDNA) has been linked to poor prognosis and higher treatment resistance in multiple types of human cancer. ecDNA are mobile genetic elements lacking centromeres that are partitioned unevenly into daughter cells at mitosis. While random segregation of ecDNA contributes to gene-copy-number heterogeneity among tumour cells, how ecDNA contribute to phenotypic heterogeneity, via mRNA and proteins, and how this affects targeted therapy outcomes is still not understood. In fact, cancer cell populations with ecDNA show remarkable heterogeneity in mRNA concentrations, which is not explained by gene-copy-number heterogeneity alone. In this paper, we use a stochastic model of gene transcription to mathematically model mRNA production and decay in cancer cells with ecDNA, extending a method that has been used to study phenotypic heterogeneity in eurokaryotic and prokaryotic cells. Using our model, we explore the space of possible mRNA distributions from a given ecDNA distribution, under different model parameter regimes. Our model is able to produce realistic mRNA distributions and we use it to estimate transcriptional parameters from single-cell measurements of ecDNA and mRNA counts in cell-line data. We discover a remarkable similarity in estimated parameters across cell lines and amplification types, indicating similar processes underlying transcription in ecDNA and chromosomal amplifications. We also model genotoxic therapy and use our model to predict the mRNA composition of persisting cells, showing the usefulness of such a modelling framework to understanding the origins of treatment resistance in ecDNA-containing cancers.

Demographic stochasticity is generally expected to contribute little to population variability when census abundance is large, because random demographic events are increasingly averaged across many individuals. This expectation, however, depends on how evenly individuals contribute to reproduction. In many natural populations, reproductive success is highly skewed, so that a large census population can coexist with a much smaller effective number of breeders (Nb). Although Nb is commonly used as a genetic measure of contemporary effective size, it also summarizes variation in reproductive contribution and may therefore contain information about demographic variability. Here, we develop a female-based age-structured theoretical framework linking Nb to temporal variability in mature population abundance, with N denoting the census number of mature females and Nb the effective number of female breeders. Recruitment variability arises from non-Poisson reproductive variation and environmental stochasticity, and its transmission to mature abundance is described by a life-history coefficient, {gamma}, determined by survival, maturation, and the contribution of overlapping cohorts. When the Poisson component of recruitment variance is negligible and Nb/N is sufficiently small, the coefficient of variation in mature abundance, CV[N], can be approximated from Nb, environmental variability, and {gamma}, with little additional information from census abundance. This approximation does not require specification of the functional form of the relationship between adult abundance and recruitment. Stochastic age-structured simulations closely reproduced the analytical prediction, even when additional demographic stochasticity in survival was included. Substantial fluctuations in mature abundance occurred at census sizes on the order of 105 when Nb was small, whereas overlapping cohorts buffered recruitment fluctuations. Because Nb can be estimated genetically from samples of a single cohort, the framework provides a direct link between genetic monitoring and population dynamics. More broadly, these results show that large census abundance alone does not guarantee weak demographic variability and highlight the potential of genetic data to inform expected population fluctuations when detailed demographic information is limited.

We propose a stochastic epidemic model to describe the spread of avian influenza in a network of poultry farms. We consider three main paths of infection transmission between poultry farms: from nearby farms, mainly through airborne diffusion; from farms belonging to the same company, through shared veterinarians, forage providers and similar; from undetected small farms or wildlife. From the parameter estimation, based on a modified Expectation-Maximization algorithm, we infer that approximately 63% of the farms were infected from nearby ones, with an infection force declining with distance; of these, more than one third belonged to the same company of the estimated infector. About 20% were infected from premises belonging to the same company but farther away than the distance threshold of 2 km; the remaining ones from wildlife or unidentified sources. These estimates have been validated by a comparison with genetic data, available for a subset of the farms: the genetic distance between two farms identified, with high probability, as an infector-infectee pair is much lower than between random pairs. These results may help in the implementation of tailored prevention and control measures for future outbreaks.

Evolution in fluctuating environments depends critically on epistatic interactions between mutations and the pleiotropic effects these mutations have on fitness across multiple conditions. However, while these effects of epistasis and pleiotropy have been extensively studied individually, the interactions between them are less well understood. In this study, we analyze how genetic background shapes the effects of mutations across multiple environmental conditions in generic models of complex traits. We show that when epistasis is widespread and random, systematic patterns of cross-environment pleiotropy can emerge naturally. Specifically we outline two regimes, widespread random epistasis (WRE) and structured pleiotropy (SP), in which patterns of global epistasis in one environment can be predicted from those in other conditions. We outline these predictions in several classes of generic null models, and test these predictions in laboratory budding yeast by measuring the fitness effects of a large set of mutations across many genetic backgrounds in several growth environments. We find that only a subset of mutations fall within the WRE and SP regimes, but that in these cases two of our models provide accurate fits to the data. Finally, we use this framework to show how modular, pleiotropic architectures shape long-term evolutionary steady states in fluctuating environments, revealing how tradeoffs between environments can develop from the geometry of shared latent interaction structure.

A major transition in the evolution of eukaryotic cells was phagocytosis: the ability of one cell to engulf another. In some cases, the engulfed cell persists, forming an endosymbiosis which is often mutually beneficial. However, it is unclear whether the interaction between the ancestors was also cooperative. There are clear benefits to developing the ability to eat a competitor, but can it ever be beneficial to eat a cooperator? We use a dynamical systems framework to identify the ecological conditions under which phagocytosis can evolve from a system of cooperating microbes.Beginning with a system of two cooperating species, we use invasion analysis to identify the parameter regimes under which a phagocytosing mutant can invade and the subsequent population dynamics. We identify a critical threshold that determines whether an invading predator will coexist with the cooperative community or replace members. Repeatedly invading with different predators, we find that the system saturates when there are three coexisting species, though this community does not have a single optimal composition. However, if we incorporate additional costs, we find there is a single evolutionarily stable community of two cooperators and a pure predator. Ultimately, we show that it is possible for phagocytosis to evolve from a cooperative community and, moreover, coexist with the original cooperating species.

Social genetic effects can substantially influence phenotypes in animals and plants and the study of them is of interest. Group housing in animals for example, in pigs provide an ideal framework for studying interactions among individuals. In animals, the genetic of social effects has been studied by simplifying the underlying interactions for example, by assuming that all the individuals in the group contribute equally to the social environment. This homogeneous behaviour is not observed in real conditions. Here, we address this issue showing that identification of interacting individuals is feasible in pig behaviour during growth, and we propose a model that leads to more accurate modelling of social genetics. Our approach infers the strength of interaction via time interval between visits of animals to automatic feeders. This measure of strength was then used in a reaction norm model to evaluate the genetic determinism of social effects. We also build a framework to transform variance components estimated for a trait continuously measured, such as feed intake per visit, into variance of total or aggregated trait e.g. during the growth period. Including social effects improves the fitting of the models. Visit-based social interaction was highest at short time intervals, decreasing rapidly and becoming negligible after 80 seconds. Variance components for the aggregated trait the total feed intake support the contribution of the social interactions. Our approach allows more precise understanding of social effects through individual feeding events and a measure of interaction strength, and further translates these into in an aggregated phenotype of interest.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Ancestral recombination graphs (ARGs) represent the coalescent history at every locus in the genome, subject to mutation and recombination. Modelling viruses using the coalescent with recombination remains less explored, due to the potential for complex, multi-scale dynamics. Using perturbative methods, we establish mathematical bounds that describe the balance between inter-host and within-host coalescence, and demonstrate that the highly-recombinant Epstein-Barr Virus (EBV) coalesces predominantly at the inter-host level. We then modify SINGER, a tool from diploid genetics, for haploid ARG inference with fine-grained variable recombination rates. We provide insight into the paradox of high recombination rates and slow linkage disequilibrium decay, as ARG-based methods potentially allow us to disentangle epistatic effects from shared ancestry, and illustrate this using the coalescent topologies of EBV's host cell-entry glycoproteins. Our findings demonstrate that ARG-based methods recover gene-level selective structure, and provide a general method for analysing virus phylodynamics in the limit of small quantities.

Bacterial mutation rates are traditionally inferred from phenotypic data using fluctuation assays. Mutation rate computation from these assays relies on a mathematical model describing the emergence of mutants during population growth. However, this standard model relies on restrictive assumptions that are often transgressed in real biological setups. Key assumptions include equal growth rate of the wild-type and the mutant population, no cell death during the growth phase, and no post-growth sampling. Although several refined mathematical models have been proposed over the last decades to circumvent some of these assumptions, none can fully account for the complex ecological scenarios frequently observed in the lab or in the wild, such as the combination of cell death and an unknown fitness effect of the mutation, or non-constant death rates. In this work, we propose to replace these standard mathematical models by a stochastic simulation framework, which can circumvent all these assumptions. We provide a fast and accurate implementation of this stochastic model, and combine it with an Approximate Bayesian Computation scheme to perform simulation-based, likelihood-free parameter inference. Our approach can simultaneously infer mutation rate and fitness effect of the mutation while accounting for death and partial plating, and could be extended to arbitrary complex conditions that can be efficiently simulated.

Species at the highest trophic levels are generally the least resilient to perturbations, yet the factors that determine their persistence remain poorly understood. We ask how the intrinsic self-regulation of a basal host governs the persistence of the trophic levels above it. To address this question, we use a minimal discrete-time model of a host-parasitoid-hyperparasitoid system motivated by the biology of trigonalid wasps, obligate hyperparasitoids that complete development only when their caterpillar host is also attacked by a primary parasitoid. We compare two forms of host self-regulation in the absence of natural enemies: scramble (Ricker) fecundity, appropriate when the parasitoid is the host's principal enemy and which recovers the classical model of Beddington and Hammond as a special case; and contest (Beverton-Holt) fecundity, appropriate when the host is held in check by many additional enemies. Using numerical continuation, Lyapunov exponents, and direct computation of the model's trajectories, we find that the form of host regulation is decisive for persistence at the top trophic level. Under scramble host dynamics, coexistence is extremely fragile with respect to both parameter values and initial conditions. Under contest host dynamics, coexistence is robust across broad parameter ranges and loses stability only gradually, through a supercritical Neimark--Sacker bifurcation and a subsequent torus-doubling route to chaos. Therefore, in at least some trophic networks, the resilience of the highest levels can be determined at the base of the network.